data=xlsread("random_coordinates.xlsx");
% 城市坐标
cities=data;
% 城市数量
nCities = length(cities);
% 参数
numParticles=50; %最大粒子数
maxIter=50;     %迭代次数
w = 1;          % 惯性权重
c1 = 2;         % 个体学习因子
c2 = 1.5;         % 全局学习因子
numVars =nCities;    % 变量的数量
T=100;  %初始温度
tem=0.99;  %温度变化率
startcity=5;     %起始城市
[bestx,bestf,bd]=sapso(cities,maxIter,w,c1,c2,numParticles,numVars,startcity,T,tem);
disp("Best Solvution");
disp(bestx);
disp("Best Function");
disp(bestf);
tusi(bd,nCities,cities,bestx);
%% 作图

function ktu=tusi(l,nCities,cities,bP)
l=[l(2:end)];
n=length(l);
count=1:n;
x=zeros(nCities+1,1);
y=zeros(nCities+1,1);
for i=1:nCities
    x(i)=[cities(bP(i),1)];
    y(i)=[cities(bP(i),2)];
end
x(i+1)=cities(bP(1),1);
y(i+1)=cities(bP(1),2);

figure(1)
plot(count,l);
xlabel('迭代次数');
ylabel('最优值');
title('模拟迭代图');
hold on
figure(2)
plot(x,y,'-ro', 'LineWidth', 2, 'MarkerEdgeColor', 'k', 'MarkerFaceColor', 'r', 'MarkerSize', 10);
xlabel('X 坐标');
ylabel('Y 坐标');
title('路径图');
grid on; % 打开网格
end


%% APSOSA
function [bestx,bestf,bd]=sapso(cities,maxIter,w,c1,c2,numParticles,numVars,startcity,T,tem)
% 距离矩阵
dist = zeros(numVars, numVars);
for i = 1:numVars
    for j = 1:numVars
        if i ~= j
            dist(i, j) = norm(cities(i, :) - cities(j, :));
        else
            dist(i, j) = inf; % 自己到自己的距离设为无穷大
        end
    end
end
city=1:numVars;%城市数
city(startcity)=[]; % 除去起始城市
% 初始化粒子位置和速度
xP = zeros(numParticles, numVars-1); % 初始化位置
xV = zeros(numParticles, numVars-1); % 初始速度为零
for i=1:numParticles
    xP(i,:)=city(randperm(numVars-1));
    xV(i,:)=randn(1,numVars-1);
    ps(i)=pathDistance([startcity,xP(i,:)],dist);%适应度
end
%个体最优(每个)
ps=ps';%适应度
pb=xP;%路径
%群体最优
[bestf,index]=min(ps);
bestx=pb(index,:);
t0=T;
wmax=w+0.5;
wmin=w-0.5;
tc1=c1;
% 迭代更新
for i=1:maxIter
   T=t0;%重置温度
   for j=1:numParticles
        % 动态调整惯性权重
        iw=wmin+(wmax-wmin)*ps(j)/mean(ps);
        r1=rand;r2=rand;
        % PSO 更新(速度)
        xV(j,:)=iw*xV(j,:)+tc1*r1.*((pb(j,:)+bestx)/2-xP(j,:))+c2*r2.*((pb(j,:)-bestx)/2-xP(j,:));
        %路径更新
        % 使用速度更新路径
        [~, sortedIndices] = sort(xV(j, :));
        newPath =xP(j, sortedIndices);
        
        % 局部搜索：2-opt算法优化路径
       newPath = localSearch2Opt([startcity,newPath], dist);
        % 去掉起始点，保持与速度一致
        newPath(1) = [];


        %更新个体最优、群体最优
        cs=pathDistance([startcity,xP(j,:)],dist);% 评估新路径的适应度
        if cs>ps(j)
            accept=exp(-(cs-ps(j))/T);%退火概率公式
            if rand()<accept
                pb(j,:)=newPath;
                ps(j)=cs;
                xP(j,:)=pb(j,:);
            end
        else
            pb(j,:)=newPath;
            ps(j)=cs;
            xP(j,:)=pb(j,:);
        end
        T=T*tem;
        if exp(-(cs-ps(j)))>rand
            tc1=c1*exp(-t0/T)*rand;
        end
        ct=pathDistance([startcity,pb(j,:)],dist);
        if ct<bestf
            bestf=ct;
            bestx=pb(j,:);
        end
        
   end
   bd(i)=bestf;
end
% 返回最优解和对应的适应度值
bestx=[startcity, bestx,startcity];
end

% 计算路径的总距离
function distance = pathDistance(path, distanceMatrix)
    n = length(path);
    distance = 0;
    for i = 1:n-1
        distance = distance + distanceMatrix(path(i), path(i+1));
    end
    distance = distance + distanceMatrix(path(n), path(1)); % 返回起点
end
% 2-opt局部搜索算法
function newPath = localSearch2Opt(path, distanceMatrix)
    n = length(path);
    newPath = path;
    improve = true;
    while improve
        improve = false;
        for i = 2:n-2 % 确保不改变起始点
            for j = i+1:n-1
                newPath = twoOptSwap(path, i, j);
                if pathDistance(newPath, distanceMatrix) < pathDistance(path, distanceMatrix)
                    path = newPath;
                    improve = true;
                end
            end
        end
    end
end
% 2-opt交换操作
function newPath = twoOptSwap(path, i, k)
    newPath = path;
    newPath(i:k) = path(k:-1:i);
end